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METHOD OF UPSCALTNG PERMEABILITY FOR 



UNSTRUCTURED GRIDS 



REFERENCE TO RELATED APPLICATIONS 



5 



This application claims the benefit of U. S. Provisional Application No. 
60/140,700 filed June 24, 1999. 



FIELD OF THE INVENTION 



10 



This invention relates generally to simulating fluid flow in a porous medium 



and, more specifically, to a flow-based method of scaling-up permeability associated 
with a fine-grid system representative of the porous medium to permeability 
associated with a coarse-grid system also representative of the porous medium. 



Numerical simulation is widely used in industrial fields as a method of 
simulating a physical system by using a computer. In most cases, there is a desire to 
model the transport processes occurring in the physical systems. What is being 

20 transported is typically mass, energy, momentum, or some combination thereof. By 
using numerical simulation, it is possible to reproduce and observe a physical 
phenomenon and to determine design parameters without actual laboratory 
experiments or field tests. It can be expected therefore that design time and cost can 
be reduced considerably. 

25 One type of simulation of great interest is a process of inferring the behavior 

of a real hydrocarbon-bearing reservoir from the performance of a numerical model 
of that reservoir. The objective of reservoir simulation is to understand the complex 
chemical, physical, and fluid flow processes occurring in the reservoir sufficiently 
well to predict future behavior of the reservoir to maximize hydrocarbon recovery. 



15 



BACKGROUND OF THE INVENTION 
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Reservoir simulation often refers to the hydrodynamics of flow within a reservoir, but 
in a larger sense reservoir simulation can also refer to the total petroleum system 
which includes the reservoir, injection wells, production wells, surface flowlines, and 
surface processing facilities. 
5 The principle of numerical simulation is to numerically solve equations 

describing a physical phenomenon by a computer. Such equations are generally 
ordinary differential equations and partial differential equations. These equations are 
typically solved using numerical methods such as the finite difference method, the 
finite element method, and the finite volume method among others. In each of these 

10 methods, the physical system to be modeled is divided into smaller cells (a set of 

which is called a grid or mesh), and the state variables continuously changing in each 
cell are represented by sets of values for each cell. An original differential equation is 
replaced by a set of algebraic equations to express the fundamental principles of 
conservation of mass, energy, and/or momentum within each smaller unit or cells and 

15 of mass, energy, and/or momentum transfer between cells. These equations can 
number in the millions. Such replacement of continuously changing values by a 
finite number of values for each cell is called "discretization". In order to analyze a 
phenomenon changing in time, it is necessary to calculate physical quantities at 

~* discrete intervals of time called timesteps, irrespective of the continuously changing 

20 conditions as a function of time. Time-dependent modeling of the transport processes 
proceeds in a sequence of timesteps. 

In a typical simulation of a reservoir, solution of the primary unknowns, 
typically pressure and phase saturation or composition, are sought at specific points in 
the domain of interest. Such points are called "gridnodes" or more commonly 

25 "nodes." Cells are constructed around such nodes, and a grid is defined as a group of 
such cells. The properties such as porosity and permeability are assumed to be 
constant inside a cell. Other variables such as pressure and phase saturation are 
specified at the nodes. A link between two nodes is called a "connection." Fluid flow 
between two nodes is typically modeled as flow along the connection between them. 
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In conventional reservoir simulation, most grid systems are structured. That 
is, the cells have similar shape and the same number of sides or faces. Most 
commonly used structured grids are Cartesian or radial in which each cell has four 
sides in two dimensions or six faces in three dimensions. While structured grids are 
5 easy to use, they lack flexibility in adapting to changes in reservoir and well geometry 
and often can not effectively handle the spatial variation of physical properties of 
rock and fluids in the reservoir. Flexible grids have been proposed for use in such 
situations where structured grids are not as effective. A grid is called flexible or 
unstructured when it is made up of polygons (polyhedra in three dimensions) having 
10 shapes, sizes, and number of sides or faces that can vary from place to place. 
Unstructured grids can conform to complex reservoir features more easily than 
structured grids and for this reason unstructured grids have been proposed for use in 
reservoir modeling. 

One type of flexible grid that can be used in reservoir modeling is the Voronoi 

1 5 grid. A Voronoi cell is defined as the region of space that is closer to its node than to 
any other node, and a Voronoi grid is made of such cells. Each cell is associated with 
a node and a series of neighboring cells. The Voronoi grid is locally orthogonal in a 
geometrical sense; that is, the cell boundaries are normal to lines joining the nodes on 
the two sides of each boundary. For this reason, Voronoi grids are also called 

20 perpendicular bisection (PEBI) grids. A rectangular grid block (Cartesian grid) is a 
special case of the Voronoi grid. The PEBI grid has the flexibility to represent 
widely varying reservoir geometry, because the location of nodes can be chosen 
freely. PEBI grids are generated by assigning node locations in a given domain and 
then generating cell boundaries in a way such that each cell contains all the points 

25 that are closer to its node location than to any other node location. Since the inter- 
node connections in a PEBI grid are perpendicularly bisected by the cell boundaries, 
this simplifies the solution of flow equations significantly. For a more detailed 
description of PEBI grid generation, see Palagi, C.L. and Aziz, K.: "Use of Voronoi 
Grid in Reservoir Simulation," paper SPE 22889 presented at the 66th Annual 

30 Technical Conference and Exhibition, Dallas, TX (Oct. 6-9, 1991). 
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The mesh formed by connecting adjacent nodes of PEBI cells is commonly 
called a Delaunay mesh if formed by triangles only. In a two-dimensional Delaunay 
mesh, the reservoir is divided into triangles with the gridnodes at the vertices of the 
triangles such that the triangles fill the reservoir. Such triangulation is Delaunay 
5 when a circle passing through the vertices of a triangle (the circumcenter) does not 
contain any other node inside it. In three-dimensions, the reservoir region is 
decomposed into tetrahedra such that the reservoir volume is completely filled. Such 
a triangulation is a Delaunay mesh when a sphere passing through the vertices of the 
tetrahedron (the circumsphere) does not contain any other node. Delaunay 
10 triangulation techniques are well known; see for example U.S. patent 5,886,702 to 
Migdal et al. 

Through advanced reservoir characterization techniques, it is common to 
model the geologic structure and stratigraphy of a reservoir with millions of grid 
cells, each populated with a reservoir property that includes, but is not limited to, rock 

15 type, porosity, permeability, initial interstitial fluid saturation, and relative 

permeability and capillary pressure functions. However, reservoir simulations are 
typically performed with far fewer grid cells. The direct use of fine-grid models for 
reservoir simulation is not generally feasible because their fine level of detail places 
prohibitive demands on computational resources. Therefore, a method is needed to 

20 transform or to scale up the fine-grid geologic reservoir model to a coarse-grid 

simulation model while preserving, as much as possible, the fluid flow characteristics 
of the fine-grid model. 

One key fluid flow property for reservoir simulation is permeability. 
Permeability is the ability of a rock to transmit fluids through interconnected pores in 

25 the rock. It can vary substantially within a hydrocarbon-bearing reservoir. Typically, 
permeabilities are generated for fine-scale models (geologic models) using data from 
well core samples. For simulation cells, the heterogeneities of the geologic model are 
accounted for by determining an effective permeability. An effective permeability of 
a heterogeneous medium is typically defined as the permeability of an equivalent 

30 homogeneous medium that, for the same boundary conditions, would give the same 
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flux (amount of fluid flow across a given area per unit time). Determining an 
effective permeability, commonly called permeability upscaling, is not 
straightforward. The main difficulty lies in the interdependent influences of 
permeability heterogeneities in the reservoir and the applied boundary conditions. 
5 Many different upscaling techniques have been proposed. Most of these 

techniques can be characterized as (1) direct methods or (2) flow-based methods. 
Examples of direct methods are simple averaging of various kinds (e.g., arithmetic, 
geometric and harmonic averaging) and successive renormalization. The flow-based 
techniques involve the solution of flow equations and account for spatial distribution 

10 of permeability. In general, the flow-based methods require more computational 
effort but are more accurate than the direct methods. 

An overview of different upscaling techniques is provided in the following 
papers: Wen, X. H. and Gomez-Hernandez, J. J., "Upscaling Hydraulic 
Conductivities in Heterogeneous Media: An Overview," Journal of Hydrology, Vol. 

15 183 (1996) 9-32; Begg, S.H.; Carter, R.R. and Dranfield, P., "Assigning Effective 
Values to Simulator Gridblock Parameters for Heterogeneous Reservoirs," SPE 
Reservoir Engineering (Nov. 1989) 455-465; Durlofsky, L. J., Behrens, R. A., Jones, 
R. C, and Bernath, A., "Scale Up of Heterogeneous Three Dimensional Reservoir 
Descriptions," Paper SPE 30709 presented at the Annual Technical Conference and 

20 Exhibition, Dallas, TX (Oct. 22-25, 1995); and Li, D., Cullick, A., Lake, L. W., 
"Global Scale-up of Reservoir Model Permeability with Local Grid Refinement", 
Journal of Petroleum Science and Engineering, Vol. 14 (1995) 1-13. The upscaling 
techniques proposed in the past were primarily focused on structured grids. A need 
exists for a method of upscaling permeabilities associated with a fine-scale geologic 

25 model to permeabilities associated with an unstructured, coarse-scale reservoir 
simulation model. 

SUMMARY 

A method is provided for scaling up permeabilities associated with a fine-scale 
grid of cells representative of a porous medium to permeabilities associated with an 
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unstructured coarse-scale grid of cells representative of the porous medium. The first 
step is to generate an areally unstructured, Voronoi, computational grid using the 
coarse-scale grid as the genesis of the computational grid. The cells of the 
computational grid are smaller than the cells of the coarse-scale grid and each cell of 
5 the computational grid and the coarse-scale grid has a node. The computational grid 
is then populated with permeabilities associated with the fine-scale grid. Flow 
equations, preferably single-phase, steady-state pressure equations, are developed for 
the computational grid, the flow equations are solved, and inter-node fluxes and 
pressure gradients are then computed for the computational grid. These inter-node 

10 fluxes and pressure gradients are used to calculate inter-node average fluxes and 
average pressure gradients associated with the coarse-scale grid. The inter-node 
average fluxes and average pressure gradients associated with the coarse grid are then 
used to calculate upscaled permeabilities associated with the coarse-scale grid. 

In a preferred embodiment, the computational grid is constructed from the 

1 5 coarse-scale grid to produce inter-node connections of the computational grid that are 
parallel to the inter-node connections of the coarse-scale grid. The cells of the 
computational grid are preferably approximately the same size as the fine-scale cells. 
The computational grid is preferably populated with permeabilities by assigning to a 
given node of the computational grid, a predetermined permeability of a cell of the 

20 fine-scale grid that would contain the given node's location if the computational grid 
were superimposed on the fine-scale grid. The flow equations that are developed for 
the computational grid are preferably single-phase, steady-state equations. The inter- 
node average fluxes and average pressure gradients associated with the coarse-scale 
grid are preferably calculated using only the inter-node connections of the 

25 computational grid that fall within a predetermined sub-domain of the computational 
grid, and more preferably such calculations are made using only the inter-node 
connections of the computational grid that are parallel to the inter-node connection of 
the coarse-scale grid associated with the sub-domain. The permeabilities associated 
with inter-node connections of the coarse-scale grid are preferably determined by 
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calculating the ratio of the inter-node average fluxes to the inter-node average 
pressure gradients that were calculated for the coarse-scale grid. 



BRIEF DESCRIPTION OF THE DRAWINGS 

The present invention and its advantages will be better understood by referring 
5 to the following detailed description and the following drawings in which like 
elements are given like numerals and letters. 

Fig. 1 illustrates a two-dimensional, structured, fine-scale grid commonly 
used in geologic modeling of a porous medium. 

Fig. 2 illustrates a two-dimensional, unstructured, coarse-scale grid having 
10 36 cells of which four cells are shown having nodes designated I, J, K and L. 

Fig. 3 illustrates a two-dimensional, unstructured computational grid 
generated from the coarse-scale grid of Fig. 2; the cells of the computational grid 
were generated using a grid-refinement ratio of 2. 

Fig. 4 illustrates 10 unstructured cells generated in the practice of this 
15 invention from three unstructured cells of Fig. 2 having nodes I, J, and K; such 10 
cells were generated from the three cells using a grid-refinement ratio of 3. 

Fig. 5 illustrates 15 unstructured cells generated in the practice of this 
invention from the three unstructured cells of Fig. 2 having nodes I, J, and K; such 
15 cells were generated from the three cells using a grid-refinement ratio of 4. 
20 Fig. 6 illustrates a two-dimensional computational grid generated from the 

grid of Fig. 2 using a grid-refinement ratio of 4, and shows a diamond shaped sub- 
domain KILL 

Fig. 7 illustrates four coarse-grid cells of Fig. 2 having nodes I, J, K, and L 
and the diamond shaped sub-domain R of Fig. 6. 
25 Fig. 8 illustrates a computational grid that is a sub-domain of the grid in Fig. 

6 and only includes the region of integration R of Fig. 6. 
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Fig. 9 illustrates a computational grid that is a sub-domain of the grid in 
Fig. 6, which includes a region of integration R of Fig. 6 and a buffer zone (skin 
region) around region R . 

The drawings illustrate specific embodiments of practicing the method of this 
5 invention. The drawings are not intended to exclude from the scope of the invention 
other embodiments that are the result of normal and expected modifications of the 
specific embodiments. 

DETAILED DESCRIPTION OF THE INVENTION 

This invention provides a flow-based method for upscaling permeabilities 

1 0 associated with a fine-scale, geologic grid system representative of a porous medium 
to estimated permeabilities associated with a coarse-scale, unstructured grid system 
representative of the same porous medium. The invention is particularly useful in 
upscaling permeabilities associated with a geologic model of a reservoir, typically 
represented by structured cells, to permeabilities associated with a coarse-scale, 

1 5 unstructured, Voronoi grid (also known as perpendicular-bisection or PEBI grid). 

As a first step in carrying out the scale-up method, an areally, unstructured, 
Voronoi grid (often referred to herein as a "computational grid") is generated using 
the coarse grid as its genesis. In accordance with the practice of this invention, the 
inter-node connections of both the coarse grid and computational grid are aligned 

20 (parallel) with each other. The computational grid is then populated with 

permeabilities from the fine-scale, geologic grid. The next step is to set up flow 
equations, preferably single-phase, steady-state pressure equations, for the 
unstructured computational grid, solve the flow equations, and compute inter-node 
fluxes and pressure gradients for the computational grid. These fluxes and pressure 

25 gradients are then used to calculate inter-node average fluxes and average pressure 
gradients associated with the coarse grid. Permeabilities associated with the coarse- 
grid connections are then calculated using the average fluxes and average pressure 
gradients calculated earlier for the coarse grid. 
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One embodiment of the invention will now be described with reference to the 
drawings. Although the drawings illustrate only two-dimensional (2-D) grid systems, 
this invention is not limited to 2-D grids. As discussed later herein, the invention can 
also be applied to three-dimensional (3-D) grids. 
5 Fig. 1 illustrates an example of a fine-scale, geologic grid for use in 

representing a porous medium such as an aquifer or hydrocarbon-bearing reservoir. 
While Fig. 1 shows an 8x8 Cartesian grid of 64 cells, it should be understood that a 
typical geologic grid could contain millions of cells. All of the cells in Fig. 1 are 
structured, which is typical of geologic grids, with each cell having a node at its 

10 center. In Fig. 1, and in the other drawings, the thick dots denote cell nodes; the 

continuous lines denote the boundaries of cells; and in Figs. 1-5 and 7 the dashed lines 
between nodes are referred to as inter-node connections. The triangles formed by the 
dashed lines in these figures form Delaunay meshes. 

Fig. 2 illustrates an example of an unstructured, PEBI grid system containing 

15 coarse-scale cells suitable for use in reservoir simulation. The grids of Figs. 1 and 2 
represent the same domain of a porous medium. In most reservoir simulation 
applications, the average size of the coarse-scale cells of Fig. 2 would be several times 
larger than the average size of fine-scale cells of Fig. 1 . A reservoir simulation model 
may for example involve from about 10 to 1000 geologic cells for each reservoir 

20 simulation cell, with the length scale changing from 10 meters or more in the geologic 
model to 100 meters or more in the reservoir simulation model. Although the cells 
depicted in Fig. 1 are smaller than the cells depicted in Fig. 2 (for sake of clarity in 
presentation) they could be even smaller in actual cases. Techniques for generating 
both fine-scale, geologic grids and coarse-scale simulation grids are well known. The 

25 inventors have discovered an effective flow-based upscaling method for transforming 
the permeabilities associated with the fine-scale grid (such as a structured grid of 
Fig. 1) to an unstructured, coarse-grid system (such as an unstructured grid of Fig. 2). 

Fig. 3 illustrates one example of a computational grid generated in the practice 
of this invention that can be used in upscaling the permeabilities from a geologic 

30 model of Fig. 1 to an unstructured reservoir simulation model of Fig. 2. The first step 
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in the method of this invention is to generate a computational grid for the same 
domain covered by Figs. 1 and 2. 

The computational grid is generated in such a way as to produce unstructured, 
PEBI cells that are smaller than the coarse-grid cells. Each inter-node connection of 
5 the computational grid is parallel to an inter-node connection of the coarse grid. As 
will be discussed in more detail below with respect to the preferred embodiment, this 
parallel alignment facilitates the upscaling of pressure gradients and fluxes which is 
an important step in the method of this invention. 

The computational grid is generated by assigning nodes in addition to nodes of 

10 the coarse-grid and regenerating PEBI cells. The additional nodes are assigned 

proportionally on all coarse-grid connections and within each Delaunay triangle. The 
same number of new nodes is added to each connection and any number of new nodes 
may be added. Increasing the number of nodes increases the number of smaller PEBI 
cells that can be formed from the coarse-grid cells. Increasing the number of smaller 

1 5 cells of the computational grid is referred to herein as increasing the refinement. The 
degree of refinement can be measured by counting the number of additional nodes per 
coarse-grid connection. If one node is added per coarse-grid connection, the 
refinement ratio is 2 since it divides the coarse-grid connection into two equal 
minigrid connections; if 2 nodes are added per coarse-grid connection, the refinement 

20 ratio is 3, and if 3 nodes are added, the refinement ratio is 4, and so on. 

The PEBI grid of Fig. 3 was generated with a refinement ratio of 2. Referring 
to Fig. 3, nodes I, J, and K represent nodes of the three cells illustrated in Fig. 2 
having nodes I, J, and K. To generate the grid of Fig. 3, one additional node was 
introduced midway between nodes K and I, one additional node was introduced 

25 midway between nodes K and J, and one additional node was introduced midway 
between nodes I and J, and so on for all other inter-node connections of Fig. 2. The 
PEBI grid of Fig. 3 was then generated based on the original nodes of Fig. 2 and the 
newly added nodes. 

Examples of computational grids having refinement ratios of 3 and 4 are 

30 shown in Figs. 4 and 5, respectively. Fig. 4 shows 10 PEBI cells generated from the 
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three cells of Fig. 2 having nodes I, J, and K, representing a refinement ratio of 3. To 
generate the PEBI cells of Fig. 4, two nodes 10 and 1 1 were added proportionately 
between nodes K and I; two nodes 13 and 14 were added proportionately between 
nodes K and J; two nodes 1 5 and 1 6 were added proportionately between I and J; and 
5 one node 12 was introduced within the triangle formed by node K, I, and J. The 

Delaunay mesh (the dashed lines) formed by inter-node connections of the 10 cells of 
Fig. 4 consists of 9 equal-sized Delaunay triangles having equal angles. All 9 smaller 
Delaunay triangles of Fig. 4 are similar to the Delaunay triangle IJK of Fig. 2. 

Fig. 5 shows 15 PEBI cells generated from the three cells of Fig. 2 having 

10 nodes I, J, and K, representing a refinement ratio of 4. To generate the 15 PEBI cells 
shown in Fig. 5, three nodes 20, 21, and 22 were added proportionately between 
nodes K and I of Fig. 2; three nodes 23, 24, and 25 were added proportionately 
between nodes K and J; three nodes 26, 27, and 28 were added proportionately 
between I and J; and three nodes 29, 30, and 31 were proportionately added within the 

1 5 triangle formed by node K, I, and J. Nodes 29, 30, and 3 1 lie at the intersections of 
lines drawn parallel to sides of triangle KIJ through nodes 20-28. Hence, all inter- 
node connections of the PEBI cells in Fig. 5, are parallel to one of the inter-node 
connections between nodes K and I, between nodes K and J, and nodes I and J of 
Fig. 2. For example, referring to Fig. 5, the inter-node connection between nodes 20 

20 and 23 is parallel to the inter-node connection between I and J of Fig. 2; the inter- 
node connection between nodes 23 and 24 is parallel to the inter-node connection 
between nodes K and J of Fig. 2; and the inter-node connection between nodes 29 and 
30 is parallel to the inter-node connection between nodes K and I of Fig. 2. 

For a given refinement ratio, the number of similar, smaller, Delaunay 

25 triangles of a computational grid that can be produced is n 2 times the number of 

Delaunay triangles of the coarse-scale grid, where n is a desired integer refinement 
ratio. For example, a computational grid with a refinement ratio of 4 will have 
sixteen (4 2 ) Delaunay triangles (the Delaunay grid of Fig. 5) that are similar to, and 
smaller than, the Delaunay triangle KIJ of Fig. 2. 
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The refinement ratio can be any integer greater than 1 . For most reservoir 
simulation applications, as a non-limiting example, the refinement ratio may range 
from 2 to about 10, and more typically the ratio will range from 2 to 5. The desired 
refinement ratio will depend on the relative sizes of the geologic grid and the coarse 
5 grid and the desired number of computational grid cells. The refinement ratio is 

preferably chosen to generate computational grid cells of approximately the same size 
as cells of the geologic grid for the same domain of the reservoir, and more preferably 
the refinement ratio chosen produces computational grid cells slightly smaller than the 
geologic cells. 

10 The method of this invention generates PEBI cells when all angles of each 

Delaunay-mesh triangle form coarse-grid connections that are smaller than 90° (acute 

o 

triangle) or when one angle is at most equal to 90 (right triangle). However, if any 

o 

angle of a Delaunay triangle is greater than 90 (obtuse triangle), a true perpendicular 
bisector does not exist. For such triangles, the grid-generation technique of this 

15 invention will produce connections that are not parallel to the overlying coarse-grid 
connections. If any such non-parallel connections are generated, they are preferably 
not used in upscaling permeabilities from the computational grid to the coarse-grid as 
discussed in more detail below. The coarse-grid is preferably generated such that 
there are no obtuse triangles in the Delaunay-mesh. This can be achieved by 

20 adjusting node locations in most situations. 

The minigrid-generation technique of this invention is not limited to a 
triangular connection mesh (Delaunay-mesh); it is also applicable to structured grids 
in which the connection mesh comprises rectangles. The present computational grid 
generation technique is also applicable to a connection mesh that comprises a 

25 combination of triangles and rectangles. 
Permeability Population 

Once the computational grid has been generated, the next step in the method 
of this invention is to populate permeabilities onto the computational grid. This 
populating step is carried out using predetermined permeabilities associated with the 
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fine-scale, geologic model (Fig. 1 for example). Various populating approaches can 
be used. One approach assigns permeabilities from the geologic grid to nodes of the 
computational grid as follows. For a given node on the computational grid, the 
geologic grid is searched until the cell of the geologic grid is located that contains the 
5 node of the computational grid, assuming that the computational grid is superimposed 
on the geologic grid. The predetermined permeability of the geologic cell that 
contains a given node of the computational grid is assigned as the permeability of that 
given node. Each inter-node connection permeability of the computational grid is 
then computed by averaging (preferably by harmonic averaging) the two node 

1 0 permeabilities forming the inter-node connection. Another approach to populating the 
computational grid with permeabilities is to assign to the midpoints of the 
computational grid's inter-node connections, permeabilities obtained from the 
geologic cells. For a given midpoint of the computational grid's inter-node 
connection, a cell of the geologic grid is located that contains the given mid-point 

1 5 (assuming that the computational grid is superimposed on the geologic grid) and the 
permeability of that geologic cell is assigned to the given midpoint. If the geologic 
model has a diagonal permeability tensor, the two areal permeabilities can be 
combined based on the connection direction, before population. 
Fluid Flow Equations 

20 After the computational grid has been populated with permeabilities, the next 

step in the practice of this invention is to develop flow equations for the 
computational grid, preferably single-phase, steady-state pressure equations for each 
cell of the computational grid, and using an assumed set of boundary conditions, solve 
the equations for each computational grid cell. The single-phase, steady-state 

25 pressure equation is given as: 

v£w> = o (1) 

where P is pressure and k is the permeability tensor. 
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The following description makes use of mathematical symbols, many of which 
are defined as they occur throughout the text. Additionally, for purposes of 
completeness, a symbols table containing definitions of symbols used herein is 
presented following the detailed description. 
5 Assuming, but not limited to, a scalar connection permeability, the discretized 

form of Eq. (1) for PEBI grids is: 

2>,fc-*)-« (2) 

1 0 where T is transmissibility , subscript j refers to the node of interest and subscript £ refers 
to all of its neighbors. The term transmissibility as used in this description refers to a 
measure of the capability of a given viscosity fluid to move across a cell boundary (or 
inter-node connection) under a pressure drop. More specifically, transmissibility is 
known to those skilled in the art as a measure of the ability of a fluid to flow between two 

kA 

1 5 neighboring cells within a porous medium. Transmissibility is expressed as — , where k 

As 

is the effective permeability of the porous medium, A is the area of the boundary between 
the neighboring cells, and Ay is the average or characteristic distance that the fluid must 
travel in moving between the two cells. 

One acceptable set of boundary conditions exerts a constant pressure gradient 

20 in the flow direction and assumes no flow in the transverse direction (no-transverse- 
flow boundary conditions). As shown in Fig. 6, fluid flow is assumed to be injected 
or produced only in the direction of the arrows (from left to right). Another set of 
boundary conditions suitable for the practice of this invention is linear boundary 
conditions.^ In linear boundary conditions, a constant pressure gradient is applied in 

25 the flow direction along all boundaries. This can result in flow across all boundaries 
which can help in upscaling permeabilities in isolated areas. The foregoing boundary 
conditions should not be construed as limitation of the invention; other suitable 
boundary conditions can also be used. 
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The pressure equation is solved in each of the two principal directions for 2-D 
applications (and in three principal directions for 3-D applications). Fluxes and 
pressure gradients are then computed for all inter-node connections of the 
computational grid. Persons skilled in the art of reservoir simulation would be 
5 familiar with development of suitable single-phase pressure equations for each cell of 
the computational grid, solving the pressure equations, and computing fluxes and 
pressure gradients for all inter-node connections of the computational grid. See for 
example, Verma, S., "A Control Volume Scheme for Flexible Grids in Reservoir 
Simulation," paper SPE 37999 presented at the Reservoir Simulation Symposium, 
10 Dallas TX (June 8-11, 1997). 
Upscaling 

Once the fluxes and pressure gradients computed are determined for the 
computational grid, the average fluxes and average pressure gradients associated with 
connections of the coarse grid are computed. These pressure gradients and fluxes are 
15 averaged over predetermined integration sub-domains associated with each coarse-grid 
cell, preferably sub-domains associated with each inter-node connection of the coarse 
grid. The ratio of the upscaled flux to upscaled pressure gradient then gives the upscaled 
permeability. This upscaled permeability can then be used to compute transmissibility. 
A mathematical basis for upscaling permeabilities from a computational grid 
20 to the coarse grid will now be provided. This description makes use of several 

mathematical symbols, some of which are defined as they occur in this description. 
Additionally, for purposes of completeness, a symbols table containing definitions of 
symbols used herein is presented following the detailed description.. 

Darcy's law for 1 -D single-phase flow in direction s is given by: 



25 



30 



= _kdP 

lids (3) 

Integrating Eq. (3) with respect to volume for computing an average flux over 
the volume of sub-domain R gives: 
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\dV 



(4) 



Assuming constant viscosity, the upscaled permeability (k ) for sub-domain 
R in direction s at the coarse-scale grid can now be defined as: 



Equation (5) can be expressed as a ratio of upscaled flux to upscaled pressure gradient 
by the following relationship, which was described in a paper by Rubin, Y. and 
Gomez-Hernandez, J. J., "A Stochastic Approach to the Problem of Upscaling of 
Conductivity in Disordered Media: Theory and Unconditional Numerical 
10 Simulations," Water Resources Research, Vol. 26, No. 4, pages 691-701, April, 1990: 



Approximating the integrals in Eq. (6) on the fine scale and substituting 
1 5 connection volume V f {Vj = AjAsj), the following expression for upscaled permeability 
(k t ) is obtained: 





(6) 



2> 



fi 





20 



The connection transmissibility can now be computed as: 




(8) 
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Equation (7) can be applied vising pressure gradients and fluxes computed 
from the pressure equation solution on the PEBI computational grid. An important 
step is to determine a sub-domain or region of integration associated with each 
coarse-grid connection over which the computational grid permeabilities are to be 
5 upscaled. 

Referring to Figs. 6 and 7, for upscaling areal connection permeabilities on 
PEBI grids, a suitable sub-domain is a diamond-shaped sub-domain R consisting of 
two triangles (KIJ and IJL) of Fig. 2. Other suitable sub-domain shapes can also be 
used. Equation (7) can now be used to compute upscaled permeability for each sub- 

10 domain using part or all of the computational grid's inter-node connections within the 
sub-domain using a predetermined selection criterion. In a simple but efficient 
approach, the upscaled permeability uses only the inter-node connections of the 
computational grid within sub-domain R that are parallel to the coarse-grid's inter- 
node connection representative of sub-domain R. In another embodiment, all of the 

1 5 computational grid's inter-node connections in sub-domain R could be used in the 
upscaling. 

Referring to Fig. 7, the computational grid inter-node connections in sub- 
domain R are used to obtain an upscaled permeability for the inter-node connection 
between nodes I and J. Since the pressure equation is solved in two principal 

20 directions in the areal domain, two upscaled permeabilities are obtained. If no- 
transverse-flow boundary conditions are used, the final upscaled permeability is 
selected based on the areal direction with the highest average pressure gradient for 
each coarse-grid areal connection. Other techniques could be applied to upscale the 
permeabilities using combinations of the two pressure equation solutions. For 

25 example, when linear boundary conditions are used, a linear combination of the two 
pressure gradients can be obtained as follows: 




(9) 
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Constants a and b are determined such that the resulting pressure gradient is 
aligned with the connection direction for which the permeability is to be upscaled. 
The upscaled permeability now can be expressed as: 



The upscaling technique of this invention can also be applied to 3-D structured 
as well as 3-D layered PEBI grids (also referred to by some as 2!/2-D PEBI grids). 
Extension to 3-D unstructured grids is also possible. The layered PEBI grids are 

10 unstructured areally and structured (layered) vertically. One way to create such grids 
is to project 2-D areal PEBI grids on geologic sequence surfaces. Areal connection 
permeabilities can be upscaled layer-by-layer or over multiple layers. For upscaling 
vertical (inter-layer) connection permeabilities, the flow equations would be solved 
with the pressure gradient in the vertical direction. Construction of layered 3-D grids 

15 is described by (1) Heinemann, Z. E., et al., "Modeling Reservoir Geometry With 

Irregular Grids," SPE Reservoir Engineering, May, 1991 and (2) Verma, S., et al., "A 
Control Volume Scheme for Flexible Grids in Reservoir Simulation," SPE 37999, 
SPE Reservoir Simulation Symposium, Dallas, TX, June, 1997. The upscaling 
method of this invention can also be extended to full 3-D PEBI grids. 

20 A very useful aspect of this invention is that the scale-up method can be 

applied to upscale permeabilities for an arbitrary number of coarse-grid connections 
simultaneously. The fluid flow equations can be solved on a computational grid 
corresponding to a large coarse-grid domain (scale-up domain), thus permitting a 
large number of permeabilities to be upscaled at once. For example, the fluid flow 

25 equations can be solved for the entire computational grid of Fig. 6 to upscale 
permeabilities for all the coarse grid inter-node connections shown in Fig. 2. 
Alternatively, the fluid flow equations can be solved on a smaller computational grid 



5 




(10) 
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corresponding to one or a few coarse-grid connections. An example of a 
computational grid for upscaling one coarse-grid connection permeability is shown in 
Fig. 8. Here, the computational-grid coincides with the diamond-shaped region of 
integration R used in upscaling permeability for coarse-grid connection IJ. 
5 The flow equations can also be solved on a larger computational grid than the 

scale-up domain thus allowing a buffer or skin around the scale-up domain whether it 
consists of single or multiple connections. An example of this is shown in Fig. 9 
where the computational grid includes the diamond-shaped region of integration R 
(KILJ) for connection I J and a buffer zone or skin region surrounding region R. 

1 0 These options offer flexibility in selecting scale-up domain size as well as 

computational-grid size for upscaling, which could be important for scale-up accuracy 
in certain applications where actual field boundary conditions are different from the 
scale-up boundary conditions. 

Although the mathematical description of the upscaling method of this 

15 invention has been presented with a single-phase, no gravity, steady-state, scalar 
permeability formulation, the method by no means is limited to these assumptions. 
Persons skilled in the art can extend the scale-up method of this invention to include 
multiphase flow, unsteady state flow, a full permeability tensor, and gravity using 
the information disclosed in this patent. 

20 Although not shown in the drawings, in practicing the present invention, the 

properties of the domain being simulated can optionally be visualized on 2-D and 
3-D unstructured grids using suitable computerized visualization equipment. The 
visualization of the properties can be helpful in analyzing application of the 
invention. In 2-D, a connection property (such as flux, pressure gradient or 

25 connection permeability) can be represented by displaying a diamond-shaped region 
associated with each connection in color from a predefined color scale. The apexes 
of this diamond-shaped region include the two node locations forming the 
connection and the end-points of the perpendicular bisector of the connection. In 
3-D, the Delaunay triangulation can be represented by a ball and stick network. 
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Node properties (such as pressure, saturation, composition or porosity) can be 
displayed by coloring the balls using a color scale. Similarly, the connection 
properties can be displayed by coloring the sticks. The color scale can display the 
magnitude of a property. For displaying a vector property (for example, pressure 
5 gradient or flux), an arrowhead can be used to represent the flow direction. 

The principle of the invention and the best mode contemplated for applying 
that principle have been described. It will be apparent to those skilled in the art that 
various changes may be made to the embodiments described above without 
departing from the spirit and scope of this invention as defined in the following 
10 claims. It is, therefore, to be understood that this invention is not limited to the 
specific details shown and described. 



Symbols 



15 



a, b = Constants, dimensionless 

2 

A = Area, ft 

k = Permeability, md 

k = Effective Permeability of Entire Grid Domain, md 

P = Pressure, psi 

T = Transmissibility, md-ft 



20 



u = Darcy Flux, ft/day 

V = Volume, ft 3 

AP = Pressure Drop, psi 

As = Connection Length, ft 

ju = Viscosity, cp 



25 



Subscripts 



30 



c = Coarse Grid 

/ = Computational grid 

/ = Connection Index 

J,/ = Node Index 

R = Sub-domain of Integration 
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